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This  work  presents  a  physics  based  circuit  model  for  calculating  the  total  energy  dissipated  into 
neutral  species  for  nanosecond  pulsed  direct  current  (DC)  dielectric  barrier  discharge  (DBD)  plasmas. 
Based  on  experimental  observations,  it  is  assumed  that  the  nanosecond  pulsed  DBD’s  which  have 
been  proposed  for  aerodynamic  flow  control  can  be  approximated  by  two  independent  regions  of 
homogeneous  electric  field.  An  equivalent  circuit  model  is  developed  for  both  homogeneous  regions 
based  on  a  combination  of  a  resistor,  capacitors,  and  a  zener  diode.  Instead  of  fitting  the  resistance  to 
an  experimental  data  set,  a  formula  is  established  for  approximating  the  resistance  by  modeling 
plasmas  as  a  conductor  with  DC  voltage  applied  to  it.  Various  assumptions  are  then  applied  to  the 
governing  Boltzmann  equation  to  approximate  electrical  conductivity  values  for  weakly  ionized 
plasmas.  The  developed  model  is  then  validated  with  experimental  data  of  the  total  power  dissipated 
by  plasmas.  ©2013  American  Institute  of  Physics.  [http://dx.doi.org/10.1063/L4792665] 


I.  INTRODUCTION 

The  need  for  improved  control  over  aerodynamic  flow 
separation  has  increased  interest  in  the  potential  use  of 
plasma  actuators.  The  inherent  advantages  of  plasma  actuator 
flow  control  devices  include:  fast  response  time,  surface  com¬ 
pliance,  lack  of  moving  parts,  and  inexpensiveness.  However, 
it  has  been  established  that  the  actuators  which  affect  the 
flow  via  directed  momentum  transfer  are  not  effective  at 
Mach  numbers  associated  with  most  subsonic  aircraft  appli¬ 
cations.  Recently,  Roupassov  et  al }  demonstrated  that  pulsed 
plasma  actuators,  in  which  energy  imparted  to  the  flow 
appears  to  effectively  control  flow  separation,  seem  to  be 
suitable  at  Mach  numbers  (M  «  0.3)  beyond  the  capabilities 
of  the  current  plasma  induced  momentum  based  approaches. 

Given  the  fundamental  differences  between  the  novel 
pulsed  discharge  approach  and  the  more  conventional  mo¬ 
mentum  based  approaches,  there  is  a  need  to  develop  an  effec¬ 
tive  and  efficient  model  for  the  energy  delivered  to  the  flow 
by  the  plasma.  Once  calculated,  that  value  can  be  input  to  a 
computational  fluid  dynamics  solver  as  an  energy  source  term 
resulting  in  a  coupled  fluid/plasma  dynamics  model. 
Multiphysics  models  of  this  type  are  required  in  order  to  study 
detailed  flow  characteristics.  However,  detailed  numerical 
simulations  involving  plasma  kinetics  are  computationally 
prohibitive  for  a  variety  of  coupled  fluid/plasma  design  prob¬ 
lems.  To  address  this  issue,  efficient  circuit  element  models 
have  been  introduced  to  approximate  the  complex  processes 
within  plasmas.  Circuit  models  such  as  those  by  Orlov2  rely 
on  empirical  constants  which  may  not  be  applicable  to  nano¬ 
second  pulsed  discharges.  To  date,  an  approximate  model  of 
nanosecond  pulsed  plasma  actuators  has  not  been  developed. 
This  paper  deals  primarily  with  establishing  a  flexible  model 
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with  relevant  physics  that  could  be  implemented  as  an  approx¬ 
imation  for  the  energy  dissipated  within  a  plasma  for  any 
pulsed  direct  current  (DC)  dielectric  barrier  discharge  (DBD) 
configuration.  Among  the  other  goals  in  this  paper  is  to  probe 
into  the  background  processes  that  occur  within  plasmas  and 
incorporate  that  knowledge  into  the  model. 

II.  LUMPED  ELEMENT  CIRCUIT  MODEL 

One  of  the  primary  assumptions  in  creating  this  model  is 
that  nanosecond  pulsed  DBD’s  can  be  approximated  by  two 
independent  regions  of  homogeneous  electric  field.  One  such 
region,  dubbed  the  “hot  spot”  is  the  region  adjacent  to  the 
powered  electrode.  This  region  makes  up  a  small  portion  of 
the  total  discharge  area  but  was  observed  to  be  an  important 
component  of  the  plasma  discharge  and  necessary  to  obtain 
agreement  with  experimentally  measured  shock  wave  dynam¬ 
ics  by  Roupassov  et  al }  The  other  region,  dubbed  the  “tail,” 
encompasses  the  rest  of  the  plasma  discharge  and  extends  to 
the  edge  of  the  dielectric.  As  both  regions  are  independent, 
the  model  presented  in  this  paper  consists  of  a  single  network 
for  each  region  containing  a  resistor,  capacitors,  and  a  diode. 

As  shown  in  Fig.  1,  circuit  elements  that  were  used  to 
model  the  plasma  include:  an  air  capacitor  Ca,  a  dielectric 
capacitor  Cj,  a  resistor  Rf ,  and  a  zener  diode  Df.  The  air  ca¬ 
pacitor  represents  the  capacitance  between  the  dielectric  sur¬ 
face  and  the  exposed  electrode.  The  dielectric  capacitor 
represents  the  capacitance  between  the  dielectric  surface  and 
insulated  electrode  and  is  proportional  to  the  thickness  of  the 
dielectric  layer.  Thus,  the  dielectric  layer  in  the  form  of  both 
its  thickness  and  the  value  of  its  dielectric  constant  plays  an 
important  role  in  determining  the  effectiveness  of  the  plasma 
actuator.  Finally,  the  zener  diode,  introduced  by  Orlov2  is 
utilized  in  the  model  to  enforce  an  energy  threshold  value 
below  which  plasma  will  not  form. 
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FIG.  1.  Electric  circuit  model  of  a  dielectric  aerodynamic  plasma  actuator. 


Since  a  uniform  charge  distribution  along  the  top  of  the 
dielectric  is  assumed,  the  typical  asymmetric  2D  plasma  actu¬ 
ator  geometry  featured  in  Fig.  1  can  be  simplified  to  a  series 
of  homogeneous  symmetric  regions.  As  shown  in  Fig.  2,  these 
regions  include:  an  anode  sheath,  “hot  spot,”  and  “tail.”  This 
assumption  results  in  a  series  of  coupled  ID  models  that 
account  for  the  chordwise  variation  along  the  actuator.  Fig.  3 
shows  the  simplified  circuit  model  within  each  homogeneous 
region. 


Dielectric  Region 


A.  Circuit 


FIG.  3.  Region  of  homogeneous  potential,  i.e.,  “hot  spot.” 


As  displayed  in  Fig.  1,  the  lumped  element  circuit  is  a 
function  of  the  two  capacitance  values,  Ca  and  C</.  In  this 
model,  the  air  is  treated  as  both  a  conductor  to  generate  a 
physical  relationship  for  the  resistance  Rf  and  a  parallel  plate 
capacitor  to  generate  Ca .  An  advantage  of  modeling  the 
plasma  as  a  conductor  in  addition  to  a  parallel  plate  capacitor 
is  that  it  generates  a  physical  relationship  for  the  resistance, 
Rf ,  a  value  that  is  traditionally  empirically  determined.  The 
air  gap  capacitor  can  be  modeled  as3 


_  G) 

~h 


where  Aa  is  the  cross-sectional  area  of  the  air  and  ha  is  the 
approximate  height  of  the  plasma  region  of  interest.  The 
height  of  the  plasma  has  been  shown  by  Roupassov  et  al }  to 
be  approximately  independent  of  applied  voltage  for  nano¬ 
second  pulsed  DBD  actuators.  As  displayed  in  Fig.  4,  Aa  is 
the  product  of  the  span  wise  length  of  the  actuator  za  and  la, 
the  chordwise  distance  from  the  exposed  electrode  to  the  end 
of  the  dielectric  region. 

The  capacitive  element  corresponding  to  the  dielectric 
can  be  modeled  as3 
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FIG.  2.  Plasma  discharge  regions  analyzed  in  this  study. 
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where  Ad  is  the  cross-sectional  area  of  the  dielectric  capaci¬ 
tive  element  and  hd  is  the  height  of  the  dielectric  barrier 
layer.  As  displayed  in  Fig.  5,  Ad  is  the  product  of  the  span- 
wise  length  of  the  actuator  za  and  dd,  the  width  of  the  dielec¬ 
tric  region. 

Treating  the  plasma  as  a  conductor,  the  resistance  for  DC 
voltage  is  proportional  to  op,Aa,  and  ha-  It  can  be  given  as3 


Rf  = 


ha 


(3) 


Starting  from  Kirchoff’s  circuit  laws,3  the  governing  dif¬ 
ferential  equation  for  the  voltage  drop  experienced  by  the  air 
gap,  AF,  is  given  by 


dAV(t) 

dt 


dV app 

dt 


— a 

Ca  +  Cd  ) 


—  K 


AV(t) 


Rf(t)(ca  +  cdy 


(4) 


r  i  if  |£|  >  Eclit 

\  o  if  \E\  <  Ecri„ 


(5) 


where  Vapp  is  the  applied  voltage  and  k  is  the  contribution 
from  the  zener  diode.  If  the  electric  field  magnitude,  given  as 
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FIG.  4.  Sketch  of  the  capacitive  air  element. 
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|  Exposed  Electrode 

Dielectric  Capacitance  Jd  1 

Covered  Electrode  J 

FIG.  5.  Sketch  of  the  dielectric  capacitive  element. 


the  physical  evolution  in  the  number  of  particles, /(v,/,  r,  t), 
defined  such  that /(v,/,  r,  t )  dv  is  the  number  of  particles  in 
a  unit  volume  located  at  point  r,  time  t,  internal  quantum 
number  /,  and  differential  velocity  range  v  +  dv.  Using  this 
distribution  function,  the  number  of  particles  at  point  r  and 
time  t  can  be  defined  as 


|£| 


(6) 


is  greater  than  the  breakdown  electric  field,2  Ecrit ,  required  to 
ionize  air,  then  k  takes  on  a  value  of  one,  otherwise  it  is  zero 
to  signify  that  plasma  has  not  formed. 


B.  Conductivity 

To  effectively  calculate  the  resistance  governed  by  Eq. 
(3),  an  expression  must  first  be  developed  for  the  electrical 
conductivity,  ap,  of  the  plasma.  This  value  is  one  that  tradi¬ 
tionally  requires  a  numerical  approach.  To  simplify  the  prob¬ 
lem  to  a  point  where  an  analytic  formulation  can  be  used, 
numerous  simplifying  assumptions  were  used  and  are 
described  in  the  following  paragraphs. 

For  any  plasma,  the  resulting  electric  current  is  com¬ 
posed  of  two  primary  terms:  the  current  from  electrons  and 
that  from  ions.  As  the  drift  velocity  of  electrons,  we,  in  a 
non-equilibrium  plasma  is  significantly  higher  than  ions,  the 
current  density  can  be  approximated  as  the  portion  from 
electrons  if  the  number  densities,  Ne  and  A/,  are  approxi¬ 
mately  the  same.  Using  a  form  of  the  generalized  Ohm’s 
Law,  the  current  density  vector,  J  and  op ,  respectively,  can 
be  written  as 


N(r,t)  =  ^2^f(v,J,r,t)dv. 


(9) 


This  distribution  function  allows  a  mathematical 
description  to  be  developed  for  the  temporal  evolution  in  the 
number  of  particles  resulting  from  particle  collisions  within 
a  control  volume.  The  time  rate  of  change  in  the  number  of 
particles  due  to  externally  applied  fields  can  be  described  as5 

Df  =f(v  +  dv,J,r  +  dr,t  +  dt)  -f(v,J,r,t) 

Dt  dt 


A  partial  differential  equation  can  be  developed  to 
describe  Eq.  (10)  using  formulas  for  the  time  rate  of  change  of 
v  and  r  established  using  the  equation  of  motion  in  the  form  of 


dv  F 

dt  m  ’ 


(11) 


dr 

dt 


(12) 


and  the  chain  rule  of  calculus  to  obtain  the  Boltzmann  ki¬ 
netic  energy  equation  given  as 


—eNewe  =  opE, 

(7) 

e(Nefie  +  NiHi), 

(8) 

where  /uf  and  tie  represent  the  ion  and  electron  mobilities, 
respectively.  Much  like  Eq.  (7),  the  electrical  conductivity 
relation  can  be  simplified  using  the  concept  of  quasineutral¬ 
ity  which  is  defined  as  having  approximate  equal  number 
densities  for  charged  particles  of  opposite  polarity.  Thus,  as 
\xe  is  typically  three  orders  of  magnitude  larger  than  /iz,  it  is  a 
good  assumption  to  approximate  the  electrical  conductivity 
as  only  coming  from  electrons  as  long  as  Ne  is  at  least  of  the 
same  order  of  magnitude  as  A/.4  Quasineutrality  is  a  typical 
assumption  that  is  valid  as  long  as  the  plasma  being  model¬ 
ing  is  far  enough  away  from  the  powered  electrode  to  avoid 
the  boundary  layer  in  plasma  physics  called  the  sheath. 

Since  a  pulsed  DC  voltage  is  assumed,  the  activation  of 
the  external  electric  field  will  follow  the  voltage  waveform 
as  a  step  function.  Thus,  two  expressions  will  be  required  for 
( 7p ,  where  the  first  is  valid  for  the  period  when  an  external 
electric  field  is  applied,  as  shown  in  Fig.  9  from  0  to  60  ns, 
and  the  second  when  the  voltage  drop  over  the  air  gap  is 
zero.  For  the  portions  of  the  voltage  waveform  that  AV  is 
zero,  the  power  is  also  zero  according  to  Ohm’s  Law  and 
thus  the  conductivity  during  this  time  is  of  no  importance. 

To  generate  a  analytic  formula  for  the  electrical  conduc¬ 
tivity,  a  distribution  function  must  be  introduced  to  describe 


Dt  dt  dr  m  dv 


(13) 


In  order  to  approximate  the  total  derivative,  a  relaxation 
time  can  be  introduced  defined  as  the  time  taken  for  the  sys¬ 
tem  to  be  reduced  to  an  equilibrium  distribution  function.  The 
tau  approximation  can  be  given  as  t  ~  (A ao ea\v\)~l ,  where 
(Jea  is  the  collision  cross  section  between  electrons  and  atoms, 
v  is  the  average  collision  velocity,  and  Na  is  the  number  den¬ 
sity  of  atoms.5  If  the  number  of  particles  within  a  control  vol¬ 
ume  is  defined  as  a  equilibrium  distribution  function,  /0,  when 
each  particle  is  at  the  same  energy  level  as  its  neighbor,  then 
the  particle  evolution  over  time  due  to  pairwise  collisions  af¬ 
ter  an  external  force  has  been  applied  can  be  given  as5 

^  =  (14) 
Dt  t  V  ' 

As  t  only  accounts  for  collisions  between  electrons  and 
neutral  atoms,  it  is  only  accurate  in  the  event  of  weakly  ion¬ 
ized  plasma.  Plasma  actuators  considered  in  this  paper  tradi¬ 
tionally  feature  a  low  degree  of  ionization,  or  simply  the 
amount  of  air  that  is  ionized,  and  thus  can  be  treated  in  a 
weakly  ionized  limit.  In  terms  of  the  momentum-transfer 
collision  frequency  which  can  be  defined  as  the  mass  cor¬ 
rected  rate  at  which  a  particle  of  a  specific  species  collides 
with  another,  the  criteria  for  a  weakly  ionized  plasma  can  be 
given  as4 
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V ei  (15) 

This  equation  requires  that  the  collision  frequency 
between  electrons  and  ions  be  much  less  than  those  between 
electrons  and  neutrals.  Thus,  if  this  requirement  is  met,  the 
collisional  occurrences  between  electrons  and  charged  par¬ 
ticles  can  be  effectively  ignored  and  the  relaxation  time 
established  is  a  good  approximation  of  the  total  time  rate  of 
change  in  the  number  of  particles  within  a  control  volume. 

After  introducing  a  relaxation  time  to  approximate  Eq. 
(13),  an  equation  of  motion  describing  the  average  velocity 
of  electrons  can  be  established  by  multiplying  by  mev  and 
integrating  over  the  electron  velocity.  An  analytic  formula¬ 
tion  for  the  average  velocity  an  electron  experiences  due  to 
an  externally  applied  field,  we,  can  be  obtained  by  assuming 
that  the  force  term  can  be  approximated  as  the  Lorentz  force 

F  =  -eE--(vxB),  (16) 

c 

where  B  represents  the  magnetic  field  and  E  represents  the 
electric  field.  As  plasma  actuators  have  no  applied  magnetic 
field,  B  can  be  set  equal  to  zero.  Therefore,  the  equation  of 
motion  for  an  electron  describing  we  can  be  given  as 


medwe(t )  mewe(t ) 

dt  t 


=  -eE{t). 


(17) 


Solving  Eq.  (17),  a  linear  first-order  differential  equation,  an 
integral  equation  is  obtained 


we(t*)  = - exp(  -  — 


m. 


£(f)exp  -  — )  dt.  (18) 


Equation  (18)  can  be  solved  in  conjunction  with  Eq.  (7) 
to  obtain  an  expression  for  the  time  varying  conductivity  and 
in  conjunction  with  Eq.  (8)  to  obtain  an  approximation  for 
the  time  varying  electron  mobility  if  the  ion  mobility  is 


neglected.  The  electrical  conductivity  and  electron  mobility, 
respectively,  can  be  given  as 


Opif) 


Ne(t*)e2 
meE(t*) , 

e 

meE(t* )  _ 


£(^)exp 

o 

•t* 

E(t)ex  p 
o 


t-f 


t  -  r 

N*o*  is* 

a  ea 


dt , 

dt  ^ 


(19) 

(20) 


where  oea  is  a  function  of  electron  energy  and  can  be 
obtained  for  various  molecules  found  in  air  from  Phelps.6 
Numerical  values  used  in  the  model  are  included  in  the 
Appendix.  The  electron  velocity  can  be  obtained  by  assuming 
a  Maxwellian  velocity  profile.  As  a  collection  of  electrons 
within  plasma  have  a  range  of  velocities,  the  Maxwellian  ve¬ 
locity  profile  represents  the  most  probable  distribution  of 
these  velocities.  Thus,  the  distribution  of  velocities, 
/(m,  v,  w),  can  be  given  by4 


/(«,  v,  w)  =  A3  exp 


1 

2 


m(u 2  +  v2  +  w2) 


kbTe 


(21) 


A3  =Ne 


m  V 

2nkbTe) 


(22) 


Using  the  non-relativistic  definition  of  kinetic  energy,  a 
relationship  can  be  established  for  the  kinetic  energy  of  an 
electron  that  is  valid  in  the  limit  |v|  «c,  where  c  is  the  speed 
of  light.  Using  this  approximation,  the  relativistic  formula¬ 
tion  of  the  kinetic  energy  can  be  approximated  as 

Ek  =  —  me 2  ~  -mv2.  (23) 

V1  -  M  2/c2 

Averaging  Eq.  (21)  and  using  the  non-relativistic  definition  of 
kinetic  energy,  the  mean  kinetic  energy  of  electrons  becomes4 


*oo 

J  — oo 

^ me(u 2  +  v2  +  w2)exp 

1  /  2  ,  2  ,  2\ 
-me(u  +  v  +  w  ) 

du  dv  dw 

kbT 

e 

>oo 

A3  exp 

J  — OO 

^ me(u 2  +  v2  +w2) 

Z77/ 

- 

- 

1 

where  kb  is  Boltzmann’s  constant.  From  the  definition  of  ki¬ 
netic  energy,  the  relationship  between  Eav  and  |v|  with  vector 
components  (m,  v,  w)  can  be  established  and  the  average  ther¬ 
mal  velocity  of  electrons  becomes 


The  required  inputs  for  Eqs.  (19)  and  (20)  include:  Ne,  the 
number  density  of  electrons,  Na,  the  number  density  of  atoms, 
and  Te ,  the  temperature  of  the  electrons.  Among  these  values, 
Na  can  be  assumed  to  be  constant  in  time  as  the  number  den¬ 
sity  of  atoms  is  significantly  higher  than  that  of  free  electrons. 

Many  other  models  incorporate  a  constant  electron 
temperature  into  their  model.2,7  Using  experimentally 
measured  values  of  reduced  electric  field,8  E/Na  vs. 
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electron  temperature  as  detailed  in  Fig.  6,  this  model  calcu¬ 
lates  a  new  electron  temperature  at  each  time  step  by  com¬ 
paring  the  E/N  value  experienced  by  the  plasma,  produced 
from  using  Eqs.  (4)-(6),  with  Fig.  6.  It  is  assumed  in  this 
paper  that  the  effects  of  an  applied  electric  field  have  an  in¬ 
stantaneous,  on  a  time  scale  much  faster  than  10-9  s,  effect 
on  electrons. 

A  time-varying  differential  equation  that  governs  Ne  can 
be  obtained  from  the  drift-diffusion  equations.  The  electron 
continuity  equation  that  governs  Ne  can  be  given  by 

^+V-rc  =  a|re|  -Pmne,  (26) 

where  Te  is  called  the  charged  species  flux,  a  is  the 
Townsend  coefficient  of  ionization,  and  /?  is  the  recombina¬ 
tion  coefficient  between  electron  and  neutral  atoms.  By  sim¬ 
plifying  the  plasma  discharge  into  a  combination  of  two 
homogeneous  regions  plus  an  anode  sheath  and  by  invoking 
the  irrotational  property  of  electric  fields,  Eq.  (26)  can  be 
simplified  by  ignoring  any  spatial  variation  in  the  number 
density  of  electrons,  i.e.,  Te  becomes 


Te^Nepe\E\.  (27) 

Assuming  the  number  densities  of  electrons  and  ions  are 
equal,  i.e.,  in  the  quasineutral  region,  the  electron  continuity 
equation  for  air,  with  a  composition  of  80%  N2  and  20%  02, 
can  be  written  as 


dNe 

dt 


<xair\Ne^eE\  -  0.80 pN2N2e  -  0.2 /J02N2e, 


(28) 


where  ocair  and  /?,  respectively,  can  be  approximated  as  hav¬ 
ing  the  form9 


where  A  and  B  are  empirical  constants  that  have  tabulated 
values  of  15  and  365,  respectively,  for  air  at  atmospheric 
pressure.9 

C.  Discharge  development 

The  final  remaining  unknown  required  to  close  the  equa¬ 
tion  system  is  to  determine  how  the  electric  potential 
changes  over  the  horizontal  length  of  the  actuator.  To  accu¬ 
rately  model  this,  it  is  important  to  incorporate  the  wall 
effects  of  the  plasma  actuator.  In  terms  of  potential  variation, 
these  wall  effects  attract  charged  particles  of  opposite  polar¬ 
ity  and  shield  charged  particles  of  the  same  polarity. 
Therefore,  for  regions  beside  the  anode  (powered  electrode 
for  positive  pulses),  an  anode  sheath  is  developed  where  an 
attraction  of  electrons  occurs  and  a  repulsion  of  positive  ions 
occur.  For  regions  above  the  dielectric  region  on  top  of  the 
grounded  electrode,  a  cathode  sheath  is  developed  where 
positive  ions  are  collected  and  electrons  are  repelled.  Fig.  7 
illustrates  the  two  predominate  sheaths  that  are  developed  in 
asymmetric  plasma  actuators. 

It  is  important  to  consider  the  effects  of  charge  collection 
and  repulsion  in  these  regions  as  such  phenomenon  can  have 
large  effects  on  the  variation  of  the  electric  potential  over  the 
chordwise  length  and  height  of  the  plasma  discharge. 

1.  Wall  effects 

As  the  model  introduced  in  this  paper  is  interested  in 
solving  for  the  amount  of  energy  the  plasma  transfers  to  neu¬ 
tral  species,  the  relative  importance  of  both  the  anode  and 
cathode  sheath  regions  needs  to  be  established.  Energy,  Q,  is 
a  function  of  electrical  conductivity  and  electric  field 
strength  and  can  be  given  as 


=  Ap  exp  —J  ,  (29) 

Pn2  =  2.8  x  10-7  >  (30) 

/?02  =  2x  1(T700,  (31) 


FIG.  6.  Plot  of  electron  temperature  vs.  reduced  electric  field. 


Q  =  gp\E2\Vvo1.  (32) 

Equation  (32)  shows  that  the  energy  transfer  from  both 
sheath  regions  is  proportional  to  the  conductivity  of  the 
plasma  within  each  region.  As  the  cathode  sheath  region  typ¬ 
ically  has  very  few  electrons  due  to  repulsion  effects,  the 
current  within  this  region  will  be  carried  by  positively 
charged  ions.  The  mobility  of  such  ions  are  significantly  less 
than  that  of  an  electron,  so  as  a  first  order  approximation  if 
the  local  charged  species  number  density  and  electric  field 
strength  are  assumed  equal,  the  conductivity  in  the  cathode 
sheath  will  be  significantly  less  than  the  conductivity  in  a 
bulk  plasma  and  the  anode  sheath,  i.e., 


FIG.  7.  Collective  wall  effects  of  the  exposed  powered  electrode  and  virtual 
electrode  (dielectric  region).  (1)  Anode  sheath  and  (2)  cathode  sheath.12 
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&s— c  ^  Gp  ~  Gs—a-  (33) 

When  combined  with  the  fact  that  the  volume  of  both 
sheath  regions  are  orders  of  magnitude  less  than  the  bulk 
plasma,  an  order  of  magnitude  approximation  to  the  plasma 
discharge  can  be  obtained  by  ignoring  the  cathode  sheath’s 
effects  for  nanosecond  pulsed  plasma  discharges.  Although 
the  cathode  sheath  is  approximately  negligible  in  terms  of 
energy  transfer,  the  higher  electric  field  strength  and  higher 
conductivity  present  in  the  anode  sheath  contain  important 
physics  necessary  for  capturing  the  energy  transferred  by  a 
plasma  discharge.  Fig.  2  shows  the  updated  “hot  spot”  and 
“tail”  regions  with  the  anode  sheath  included  in  the  hot  spot 
region  adjacent  to  the  powered  electrode. 

2.  Electric  potential  variation 

To  establish  a  variation  in  the  electric  potential,  the  gov¬ 
erning  Maxwell  equations  can  be  used  which  are  written  as3 


V  •  E  =  — 

(34) 

eo  ’ 

V  •  B  =  0, 

(35) 

8B 

VX£  =  --, 

(36) 

„  „  dE 

V  xB  =  n0J  +  n0e 0  — , 

(37) 

in  differential  form  where  py  =  e{rii  —  ne)  is  the  net  charge 
density  and  p0  is  the  permeability  of  free  space.  In  the  ab¬ 
sence  of  a  time- varying  magnetic  field,  Eq.  (36),  Faraday’s 
Law  of  Induction,  simplifies  to 

Vxf  =  0,  (38) 

and  since  the  curl  of  E  is  zero,  the  electric  field  can  be  solved 
for  as  a  potential  function  </>  and  substituted  into  Gauss’ 
Law,  Eq.  (34).  The  resulting  equations,  respectively,  can  be 
given  by 


> 

I 

II 

(39) 

v20  =  -(X. 

(40) 

e 


The  net  charge  density  within  the  quasineutral  region  of  a 
plasma  is  equal  to  zero  as  ne  =  rii.  For  the  anode  sheath,  the 
net  charge  density  can  be  approximated  if  the  plasma  is 
assumed  to  uniformly  distribute  its  charge  density.  Therefore, 
the  quasineutral  region  of  the  “hot  spot”  and  “tail”  features 


=  < 


e  2 
XT  ■ 


eNt 
2eo 
eNP 


'  eNe  2  eNe 

.2^3°  “V1 


2eo  L 


n  2  |  ^ break  V app 

3  j 


equal  number  densities  of  both  electrons  and  ions  (ne  =  «/), 
while  the  cathode  and  anode  sheaths  only  have  electron  and 
ion  densities,  respectively.  If  this  is  assumed  then  the  anode 
sheath  can  be  approximated  as  having  a  net  charge  density 
equal  to  the  quasineutral  region’s  calculated  electron  number 
density.  By  assuming  the  anode  sheath  is  devoid  of  any  posi¬ 
tive  ions  and  has  a  time-dependent  number  density  of  elec¬ 
trons,  Poisson’s  equation  for  the  sheath  can  be  given  by 


d 2(j>  _  eNe(t) 
dx2  e0 


(41) 


If  this  form  is  assumed  for  a  single  Debye  length  {2d),  then 
the  remaining  discharge  (rest  of  “hot  spot”  and  “tail”)  region 
is  part  of  the  quasineutral  bulk  plasma  and  can  be  given  by 
Laplace’s  equation 


(42) 


The  ordinary  differential  equations  for  electric  potential  vari¬ 
ation  in  the  hot  spot  and  tail  regions,  respectively,  can  be 
solved  to  provide  an  approximate  ID  spatial  variation 


— — x2  -I-  C\X  C2  if  x  <  Xd 

2e0  (43) 

C3X  +  C 4  if  x  > 


Equation  (43)  requires  a  total  of  4  boundary  conditions. 
Those  can  be  summarized  as 


^(x  =  0)  =  Vapp,  (44) 

02  (x  =  L)  =  Vbreak,  (45) 

0j(x  =  Ad)  =  02(x  =  Ad),  (46) 

^1(x  =  Ad)=^(x  =  Ad),  (47) 

dx  dx 

where  (j)l  is  the  potential  in  the  anode  sheath,  (j)2  is  the  poten¬ 
tial  in  the  quasineutral  region,  L  is  length  of  the  actuator,  and 
2d  is  a  Debye  length.  The  first  boundary  condition  is  Vapp  at 
the  cathode  and  the  second  is  based  on  the  experimental  obser¬ 
vations  by  Roupassov  et  al. 1  It  was  observed  that  a  plasma  dis¬ 
charge  could  be  approximated  as  stopping  at  the  edge  of  the 
grounded  electrode  for  asymmetric  actuators  independent  of 
the  applied  voltage;  so,  the  edge  represents  the  absolute  limit 
of  ionization  or  the  breakdown  voltage  of  air.  Using  Eqs. 
(44)-(47)  as  boundary  conditions,  the  potential  becomes 


V break  V c 


app 


■vt 


app 


if  x  <  2d 


x 


{x  L)  H-  Vbreak  if  X  >  2d  • 


(48) 
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The  Debye  length  provides  an  order  of  magnitude 
approximation  for  the  extent  of  a  plasma  sheath  by  assuming 
an  exponential  Boltzmann  distribution  in  the  charge  density 
within  the  plasma  discharge.  Substituting  this  into  Poisson’s 
equation, 


d2(j)  ,T 

e°i?  =  eN° 


(49) 


where  is  the  charged  particle  density  far  away  from  the 
electrode.  Taking  a  first-order  Taylor  expansion,  the  Debye 
length  can  be  given  as4 


Xd  ~ 


/  ep  kbTe\2 
\Nooe2)  ’ 


(50) 


Although  at  high  voltages,  a  first-order  approximation  fails, 
Eq.  (50)  still  provides  an  order  of  magnitude  approximation 
of  the  extent  of  the  anode  sheath. 


D.  Numerical  procedure 

When  solving  for  the  energy  imparted  to  neutral  species, 
a  coupled  equation  system  results  from  Eqs.  (4)  and  (28). 


The  coupled  terms  include 

Rf(xNe,  (51) 

re  oc  AV,  (52) 

p  oc  AV,  (53) 

01  OC  Ne.  (54) 


To  solve  the  resulting  equation  system,  the  Dormand-Prince 
Runge-Kutta  method  was  employed.  This  method  provides 
an  efficient  way  to  incorporate  an  adaptive  step  size  that  is 
important  for  computational  efficiency  in  a  problem  that 
requires  small  time  steps  for  convergence.  The  benefit  of 
such  a  procedure  can  be  illustrated  through  a  simplistic 
example.  If  the  error  of  each  time  step  is  defined  as 


4  =  y(t)  -  4-  (55> 

then  if  two  step  sizes  are  considered,  h\  and  the  error  of 
each  iteration  and  their  relative  error,  respectively,  can  be 
given  as 


y(t)  -yJi  =  4i  =  ahi>  (56) 

y(t)  -  yn2  =  en2  =  ah2,  (57) 

y^-yii  =  a(hi  ~h2)-  (58) 

Therefore,  for  a  given  error  tolerance,  e,  a  sequence  of  step 
sizes  can  be  generated 


TABLE  I.  Experimental  parameters.1 


ha 

0.4  mm 

hd 

0.3  mm 

Q 

2.7 

1 

30  mm2 

30  mm2 

V 

50  kV 

AT 

60  ns 

which  allows  a  numerical  ODE  solver,  such  as  those  employ¬ 
ing  the  Dormand-Prince  method  to  minimize  functional  error 
by  adjusting  the  step  size  after  each  time  step.10  Numerical 
integration  required  for  Eqs.  (19)  and  (20)  and  energy 
derived  from  the  circuit  model,  given  as 


Q 


*  AV2 

o  R/(0 


dt , 


(60) 


were  performed  via  the  Gauss-Kronrod  quadrature  method10 
at  each  time  step. 


III.  RESULTS 

A.  Validation  against  experiment 

To  validate  the  accuracy  of  the  model  described  in  this 
paper,  comparisons  with  data  presented  in  Roupassov 
et  al. 1  are  provided.  The  experimental  parameters  that  were 
mentioned  and  used  in  the  circuit  model  are  given  in  Table 
I.  Reference  1  uses  the  electrode  configuration  detailed  in 
Fig.  8. 

Fig.  9  illustrates  an  approximation  of  the  applied  voltage 
square  wave  that  was  introduced  in  the  experimental  work  of 
Roupassov  et  al. 1  The  slope  that  is  introduced  is  to  simulate 
a  function  that  is  differentiable.  This  is  needed  for  V'app{t)  in 
the  governing  differential  equation  as  detailed  by  Eqs.  (4) 
and  (5).  One  could  also  generate  a  continuous  function  using 
Fourier  decomposition  of  a  traditional  square  wave;  how¬ 
ever,  this  would  not  account  for  the  minor  rise  and  fall  times 
found  in  experimentation. 


k+2  =  q 


{hj  ~  hj+ i)e 

\yn\X\  -y«i I’ 


FIG.  8.  Experimental  scheme  used  by  Roupassov  et  al 3  for  the  discharge 
(59)  gap.  (1)  High-voltage  electrode;  (2)  dielectric  layer;  (3)  low-voltage  elec¬ 
trode,  (4)  zone  of  discharge  propagation,  and  (5)  insulating  plane. 
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FIG.  9.  Plot  of  the  input  voltage,  Vapp  vs.  time  used  in  the  model. 
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FIG.  10.  Plot  of  electron  number  density  vs.  time  for  the  “hot  spot.” 


B.  “Hot  spot”  results 

For  the  small  region  of  0.4  mm  x  0.4  mm  over  all  time 
outside  of  the  anode  sheath,  the  time  variation  in  the  number 
density  can  be  established.  As  an  initial  condition  for  Eq. 
(28),  1015  m-3  electrons  were  assumed  based  on  work  by 
Ref.  1.  As  displayed  in  Fig.  10,  there  is  a  large  gradient  that 
occurs  during  the  rise  time  that  peaks  around  1.13  x  1019 
m-3.  It  is  also  evident  in  Fig.  10  that  the  recombination  of 
electrons  is  largely  negligible  on  the  nanosecond  time  scale. 
If  a  frequency  of  1  kHz  is  used,  the  recombination  of  elec¬ 
trons  with  atoms  allows  a  steady- state  electron  number  den¬ 
sity  to  be  achieved  on  the  nanosecond  time  scale.  The 
recombination  of  electrons  becomes  a  significant  quantity 
when  exploring  the  dynamics  of  a  plasma  discharge  on  the 
microsecond  time  scale. 

Using  Eq.  (4)  and  Fig.  6,  the  model  was  able  to  produce 
a  time-varying  electron  temperature.  As  displayed  in  Fig.  11, 
there  is  an  initial  spike  in  the  electron  temperature  to  29  eV 
(^340  000  K)  that  coincides  with  the  peak  in  electron  num¬ 
ber  density  at  1 1  ns.  Fig.  1 1  also  suggests  that  the  assumption 
of  a  constant  electron  temperature  is  not  an  accurate  assump¬ 
tion  for  nanosecond  pulsed  DBD  plasmas.  The  variability  in 
the  electron  temperature  beyond  40  ns  is  due  to  the  numerical 
error  tolerances  that  are  selected  when  solving  Eqs.  (4),  (19), 
and  (28).  Fig.  11  shows  that  when  reducing  the  relative  error 
in  the  Dormand-Prince  method  from  10-2  to  10-3,  the  strong 
functional  variability  experiences  a  significant  reduction. 
The  higher  the  gradients  are  during  the  rise  and  fall  time,  the 


finer  the  relative  error  tolerances  are  required  to  be  to  guar¬ 
antee  convergence  of  a  solution. 

Using  the  results  displayed  in  Figs.  10  and  11,  the  total 
power  dissipated  to  neutral  species  as  a  function  of  time  can 
be  established.  Using  the  result  of  Eq.  (4)  and  the  relation¬ 
ship  for  the  instantaneous  power,  the  time-varying  power 
imparted  to  the  flow  can  be  given  as7 


p(t) 


A  V2(t) 

Rfit)  ' 


(61) 


As  shown  in  Fig.  12,  the  instantaneous  power  is  domi¬ 
nate  during  the  rise  time  for  the  “hot  spot”  region.  Upon  inte¬ 
grating  the  instantaneous  power  over  time  using  Eq.  (60), 
this  model  produces  an  energy  value  of  2.1  mJ  for  this 
region.  When  compared  to  the  experimentally  determined 
value  of  4.2  mJ  by  Roupassov  et  al.}  this  model  produces  an 
order  of  magnitude  estimate  for  the  energy  imparted  to  neu¬ 
tral  species  in  this  region. 


C.  “Hot  spot”  sheath  results 

Using  the  Debye  length  approximation  provided  by  Eq. 
(50)  and  the  solution  from  the  quasineutral  “hot  spot”  region 
(ne  ~  1019  m-3,  Te  ~  10  eV),  a  Debye  length  of  7.43  /un  is 
generated.  The  electron  number  density  experienced  in  the 
anode  sheath  is  assumed  to  be  equal  to  the  values  calculated 
in  the  adjacent  “hot  spot”  region.  Fig.  13  shows  the  time 
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FIG.  11.  Plot  of  electron  temperature  vs. 
time  for  the  “hot  spot”:  (a)  10-2  error 
and  (b)  10-3  error. 
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FIG.  12.  Plot  of  power  vs.  time  for  the  “hot  spot.” 


variation  in  the  electron  number  density  within  the  sheath 
and  also  shows  a  peak  number  density  of  1.13  x  1019  m-3. 

Using  the  revised  form  of  the  electric  potential  within 
the  anode  sheath  given  by  Eq.  (48)  and  using  Fig.  6,  the  elec¬ 
tric  temperature  variation  over  the  duration  of  the  pulse  can 
be  established.  As  displayed  in  Fig.  14,  there  is  an  initial 
spike  in  the  electron  temperature  to  (^348  000  K)  that  coin¬ 
cides  with  the  peak  in  electron  number  density  at  1 1  ns  much 
like  the  hot  spot  region. 

As  shown  in  Fig.  15,  the  instantaneous  power  is  dominant 
during  the  rise  time  for  the  anode  sheath  region.  Upon  inte¬ 
grating  the  instantaneous  power  over  time  using  Eq.  (60),  this 
model  produces  an  energy  value  of  0.045  mJ  for  this  region. 
This  number  is  quite  small  compared  to  the  2.1  mJ  experi¬ 
enced  in  the  “hot  spot”  region.  However,  the  anode  sheath 
does  have  a  slightly  higher  linear  energy  density  than  the  “hot 
spot”  region  (6J/m  vs.  5.25  J/m).  The  reason  that  there  is  not 
significant  deviation  predicted  in  the  anode  sheath  and  quasi¬ 
neutral  “hot  spot”  regions  is  that  explicit  charge  buildup  is  not 
accounted  for  in  the  circuit  model  within  the  sheath  region. 
As  shown  by  Ref.  1 1,  the  electric  field  in  the  sheath  and  adja¬ 
cent  quasineutral  regions  are  approximately  equal  until  charge 
buildup  is  allowed  within  the  anode  sheath  during  the  length 
of  the  pulse  or  a  series  of  pulses. 


Time  [ns] 


FIG.  13.  Plot  of  electron  number  density  vs.  time  for  the  “hot  spot”  sheath. 


FIG.  14.  Plot  of  electron  temperature  vs.  time  for  the  “hot  spot”  sheath. 


FIG.  15.  Plot  of  power  vs.  time  for  the  “hot  spot”  sheath. 

D.  “Tail”  results 

For  the  region  0.4  mm  x  4.6  mm,  the  time  variation  in 
the  number  density  can  be  established.  As  an  initial  condi¬ 
tion  for  Eq.  (28),  1015  m-3  electrons  were  assumed,  the  same 
number  of  electrons  assumed  for  the  “hot  spot”  region. 
When  comparing  Figs.  10  and  16,  the  tail  region  experiences 
a  lower  growth  rate  in  the  number  of  electrons  which  is  due 
to  the  lower  electric  field  experienced  by  this  region.  As 
described  by  Eq.  (28),  a  lower  electric  held  produces  a  lower 


FIG.  16.  Plot  of  electron  number  density  vs.  time  for  the  “tail.” 
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FIG.  17.  Plot  of  electron  temperature  vs.  time  for  the  “tail.” 


number  of  ionizations  and  therefore  a  more  gradual  rise  and 
lower  total  peak  in  Ne ,  approximately  3.4  x  1017  m-3. 

Using  Eq.  (4)  and  Fig.  6,  the  model  was  able  to  produce 
a  time-varying  electron  temperature.  As  displayed  in  Fig.  17, 
the  traditional  assumption  of  1  eV  (^11  600  K)  does  not 
agree  well  with  the  results  obtained  in  this  model  for  the  tail 
region  on  the  nanosecond  time  scale.  Instead,  Fig.  17  sug¬ 
gests  that  the  peak  electron  temperature  is  achieved  during 
the  rise  time  of  the  pulse,  220  000  K  (^19  eV)  and  then 
trends  downward  during  the  plateau  portion  of  the  voltage 
waveform.  Much  like  Fig.  11,  the  highest  electron  tempera¬ 
tures  are  achieved  during  the  initial  high  gradient  of  the 
pulse.  Fig.  17,  unlike  Fig.  11,  also  shows  a  rise  in  electron 
temperature  during  the  fall  time  of  the  pulse  as  well.  This  is 
due  to  the  lower  electric  potential  and  negative  gradient 
experienced  in  the  tail  region  during  the  fall  time. 

As  shown  in  Fig.  18,  the  instantaneous  power  is  dominant 
during  the  plateau  portion  of  the  applied  voltage  pulse  for  the 
“tail”  region.  Upon  integrating  the  instantaneous  power  over 
time  using  Eq.  (60),  this  model  produces  an  energy  value  of 
6.6  mJ.  When  compared  to  the  experimentally  determined 
value  of  8  mJ  by  Roupassov  et  al.1  this  model  produces  an 
absolute  error  of  «17.5%.  The  significant  improvement  in 


TABLE  II.  Comparison  between  calculated  and  experimentally  measured 
energy  deposition. 


Circuit  model  (mJ) 

Experimental1  (mJ) 

Abs.  error  (%) 

Hot  spot’ 

2.1 

4.2 

50 

Tail’ 

6.6 

8 

17.5 

Total 

8.7 

12.2 

28.7 

accuracy  for  the  tail  region  can  likely  be  attributed  to  the 
larger  distance  from  the  cathode.  This  increase  in  distance 
improves  the  assumptions  of  quasineutrality  and  that  the  spa¬ 
tial  diffusion  of  charged  species  is  negligible.  The  region 
close  to  the  cathode  features  complicated  ion  and  electron 
buildup  and  as  the  relative  distance  from  the  cathode 
increases,  its  impact  on  the  problem  becomes  negligible. 

Table  II  summarizes  the  results  obtained  using  the  cir¬ 
cuit  model  presented  in  this  paper  and  associated  experimen¬ 
tal  measurements  made  by  Roupassov  et  al.1  for  both  the 
“hot  spot”  and  “tail”  regions. 

IV.  CONCLUSION 

A  new  lumped  element  circuit  model  was  presented  that 
is  valid  for  pulsed  DC  DBD  plasmas.  The  model  approxi¬ 
mates  the  total  energy  dissipated  into  neutral  species  using  a 
lumped  element  circuit  while  containing  relevant  plasma 
physics  in  the  form  of  a  variable  electron  temperature  and 
number  density.  An  approximate  expression  was  formulated 
using  the  conductivity  of  the  discharge  to  calculate  the  resist¬ 
ance  value  for  the  air  gap.  Asymmetric  wall  effects  were 
also  approximated  in  the  model  by  including  the  effect  of  the 
anode  sheath.  Results  of  the  model  were  verified  against  a 
pulsed  DC  experiment  conducted  by  Roupassov  et  al.1  and 
order  of  magnitude  agreement  was  obtained  for  the  energy 
imparted  into  the  plasma  in  both  the  homogeneous  “hot  spot” 
region  and  “tail”  region. 
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APPENDIX:  CROSS-SECTION  DATA 

The  momentum  cross  sections  used  for  the  numerical 
approximations  in  the  model  are  shown  in  Table  III.  These 
values  were  obtained  by  A.  V.  Phelps.6 


TABLE  III.  Momentum  cross  sections.6 


Kinetic  energy  (eV) 

N2  (10-16cm2) 

02  (10-16cm2) 

1 

10 

7.2 

2.1 

27 

6.6 

3 

21.7 

5.7 

4 

12.6 

5.5 

5 

10.9 

5.6 

10 

10.4 

5 

15 

11 

oo 
o 6 

20 

10.2 

8.6 

FIG.  18.  Plot  of  power  vs.  time  for  the  “tail.” 
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